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The stationary Josephson effect in a system of ballistic graphene is studied in the framework of quasiclas- 
sical Green's function theory. Reflecting the ultimate two-dimensionality of graphene, a Josephson junction 
involving a graphene sheet embodies what we call a planar Josephson junction, in which superconducting 
electrodes partially cover the two-dimensional graphene layer, achieving a planar contact with it. For cap- 
turing this feature we employ a model of tunneling Hamiltonian that also takes account of the effects of 
inhomogeneous carrier density. Within the effective mass approximation we derive a general formula for the 
Josephson current, revealing characteristic features of the superconducting proximity effect in the planar 
Josephson junction. The same type of analysis has been equally applied to mono-, bi- and arbitrary N- layer 
cases. 
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1. Introduction 

When two superconductors are placed spatially apart 
but weakly coupled, a finite amount of dissipationless 
equilibrium current is generally induced between the 
two superconductors. Though this effect, the stationary 
Josephson effect, was originally predicted 1 ' for a system 
of tunneling junction, i.e., two superconductors separated 
by a thin insulating barrier, the same effect is now known 
to exist in a broad range of superconducting junctions, 
and especially in the ones mediated by various other non- 
superconducting elements, such as a normal metal, 2 ' a 
two-dimensional electron gas, 3 ' or a quantum dot. 4 ' Yet, 
the characteristic behavior of the dissipationless current 
in such Josephson junctions is strongly influenced by 
the electronic property of the element inserted between 
the two superconductors and by the way that element is 
placed between them. 

During the last several years, a variant of such 
Josephson junctions realized in a system of graphene, 
i.e., the superconductor-graphene-superconductor (SGS) 
junction, has become a target of intense theoretical 5-11 ' 
and experimental 12-17 ' studies. Much has been studied 
on how the Josephson current is affected by the unique 
band structure of a graphene monolayer, 18 ' 19 ' in which 
the conduction and valence bands touch conically at K + 
and K_ points in the Brillouin zone (the Dirac points). 
An interesting result is that in the SGS Josephson junc- 
tions involving a monolayer of graphene, the critical cur- 
rent I c at zero temperature remains finite at e = (when 
the chemical potential is placed at the level of the Dirac 
points) in spite of the vanishing density of states. 6 ' 

On contrary, most of the theoretical studies performed 
so far have neglected the unique structural character of 
the SGS Josephson junction. Since graphene is an iso- 
lated ideal two-dimensional electron system, a natural 
way to fabricate an SGS junction is to deposit supercon- 
ducting electrodes on top of a graphene flake. 12-17 ' Then, 
the graphene sheet beneath the superconducting elec- 
trodes acquires a two-dimensional, planar contact with 



the electrodes. This is indeed a very unique situation, 
opposing, e.g., to the case of a two-dimensional electron 
gas imbedded in a semiconductor hetero-structure that 
has a one-dimensional (linear) contact with supercon- 
ducting electrodes. 3 ' In most of the existing theoretical 
studies 6 ' 7,9 ' 11 ' the planar character of the SGS junction 
is poorly taken into account. It is usually assumed that 
an energy-independent effective pair potential A e g is in- 
duced inside the graphene sheet in the region covered 
by the superconductors. But this very assumption re- 
duces the SGS junction of intrinsically planar nature to 
a conventional linear junction. The latter consists of a 
graphene sheet of a finite length placed between two (hy- 
pothetical) graphene superconductors. 

The planar character of the SGS junction manifests 
in the temperature (T-) dependence of the Joseph- 
son current. The conventional formulation based on the 
energy-independent A c fj fails, by its construction, to de- 
scribe this feature. Of course, one can always assume 
a T-dependence of A e g = A e ff(T) and discuss the In- 
dependence of Josephson current, but this is not more 
than an ad hoc solution. For example, the T-dependcncc 
of the Josephson current has been calculated by assum- 
ing A e ff(T) as described by the BCS theory. 11 ' Here, in 
this paper, extending an earlier work of ref. 20, we give 
a more fundamental solution to this issue, allowing for 
a reliable prediction of the T-dependence of Josephson 
current in the planar SGS junction. In the approach un- 
dertaken we explicitly take account of the tunneling of 
electrons in the planar regions of contact. In other words, 
we describe the coupling between the graphene sheet and 
the superconducting electrodes by a tunneling Hamilto- 
nian. The strength of the tunnel coupling is controlled 
by the parameter T. 

Another advantage of our approach, which we aim at 
reporting in this paper, is that it allows for a system- 
atic generalization to the case of multi-layer graphene. 
A renewed insight into the use of quasiclassical Green's 
function 21,22 ' in the system of SGS junction enables this. 
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Here, the electronic property of mono-, bi- and multi- 
layer graphenes is treated under the effective mass ap- 
proximation. On one hand, most of the published the- 
oretical results have focused on the case of monolayer 
graphene, with the exception of ref. 10 treating the case 
of bilayer. The experiments are, on the other hand, not 
restricted to the case of monolayer. 12 ' 13, 16 ) The bilayer 
and multi-layer (AT-layer with N — 3, 4, • • • ) cases are 
(and will be in the near future) equally highlighted exper- 
imentally. A bilayer graphene has a nearly quadratic en- 
ergy dispersion. 23 * 1 In the case of A^-layer graphene with 
N > 3, mono- and bi-layer type dispersions coexist, but a 
basic tendency is determined only by the parity of N. 24 ^ 
We show that the Joscphson effect in the SGS junction 
also exhibits such an even-odd feature with respect to N , 
the number of layers. 

The backbone of our approach is the use of a tunneling 
Hamiltonian at the level of the modeling (of planar struc- 
ture) , and of the quasiclassical Green's function approach 
on which subsequent studies of the Josephson current 
are entirely based. In addition to the full formulation of 
the present approach that has not been done before, the 
new ingredients first taken into account in this paper are 
the following. The effect of inhomogeneous carrier den- 
sity: due to the contact with superconductors, the carrier 
density in the graphene sheet becomes higher in the re- 
gion covered by superconductors. Since this results in a 
mismatch in the Fermi wave number at the interface be- 
tween the covered and uncovered regions, we expect a 
reduction of the Josephson current. This effect was ig- 
nored in ref. 20. As mentioned earlier, we also first deal 
with the case of an arbitrary number A^ of layers. Taking 
these two new elements into account, we study the T- 
and T-dependences of the Josephson current. To avoid 
unnecessary complication the chemical potential fi is as- 
sumed to be away from the Dirac point, and we consider 
only the ballistic regime. 

The paper is arranged as follows. In §2, we start by 
modeling the SGS planar Josephson junction considering 
generally the case of A^-layer graphene. We then intro- 
duce a thermal Green's function adapted for the system 
in consideration. In §3, we derive a general formula for 
the Josephson current in the cases of mono- and bi-layer 
within the quasiclassical approximation applied with the 
use of the thermal Green's function introduced in §2. 
In §4, we extend the argument in the previous section 
to the case of arbitrary AT-layer. In §5, the behavior of 
the Josephson critical current is studied in the short- 
junction limit. We discuss analytic expressions obtained 
in some specific cases of parameters. Otherwise, the crit- 
ical current is evaluated numerically. Implications of the 
obtained results are discussed in the light of current ex- 
perimental studies that focus on the same limit. Section 
6 is devoted to summary. We set fcs = h = 1 throughout 
the paper. 

2. Modeling of the SGS Junction 

We start by modeling the SGS junction, taking prop- 
erly into account both its structural uniqueness and 
the nature of electronic states in the intermediate non- 
superconducting element, i.e., in the graphene sheet. 




Fig. 1. Planar junction geometry: Josephson junction consisting 
of a graphene sheet on which two superconductors Si and S2 of 
width W are deposited with separation L. The pair potential is 
assumed to be Ae'^/ 2 in Si and Ae~ i<p / 2 in S2. 



Here, we consider the case of a graphene sheet consisting 
of generally A" layers (treating the cases of mono-, bi-, 
and multi-layer graphene in parallel), which we describe 
in the effective mass approximation. As a basic theoret- 
ical tool indispensable for the subsequent analyses, we 
then introduce a thermal Green's function adapted for 
the system of SGS junction incorporating an A^-layer 
graphene sheet. The influence of superconducting elec- 
trodes on quasiparticles in graphene is taken into account 
by a self-energy associated with that Green's function. 

The SGS junction we consider has a construction as 
depicted in Fig. 1. Two superconductors Si and S2 (of 
width W) are placed (with separation L) on top of a 
clean graphene sheet, where Si and S2 occupy the region 
of L/2 < x and that of x < —L/2, respectively. Note 
that only the top (first) layer is in contact with Si and 
S2. We assume that the pair potential in Si and S2 is 
given by 

r Ae^ 2 (L/2<x) 
A(x) = I {-L/2 <x< L/2) . (1) 
[ Ae-^/ 2 (x<-L/2) 

Let us assume that the coupling of the graphene sheet 
and the superconductors is described by a tunneling 
Hamiltonian. Then, the resulting proximity effect on 
quasiparticles in graphene is described by the self- 
energy 25 ) for the thermal Green's function as given be- 
low. 

The coupling with Si and S2 also induces carrier dop- 
ing in the graphene sheet, i.e., the carrier density in the 
covered region of |x| > L/2 becomes higher than that in 
the uncovered region of \x\ < L/2. This effect cannot be 
described by the self-energy, so we take this into account 
by adding the effective potential of a negative value — U 
only in the covered region but for every layer. This po- 
tential rapidly varies across the interface at x = ±L/2 
with a characteristic length scale which is comparable to, 
or shorter than, the Fermi wave length but longer than 
the lattice constant a of graphene. 26 ) It is convenient to 
introduce the renormalized chemical potential /x defined 
by 

fj, (-L/2<x<L/2) 
fi + U (L/2<\x\) 



(2) 
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Let us turn our attention to the electronic property 
of the graphene sheet, which we assume to be composed 
of N layers. The unit cell of the hexagonal lattice of 
monolayer graphene contains A and B sites. We assume 
that graphene layers are subjected to the AB stacking. 
That is, A sites in the first layer are located just above 
B sites in the second layer, and B sites in the second 
layer are located just above A sites in the third layer, 
and so forth. To describe electron states in an A-layer 
graphene sheet, we first employ a tight-binding model 
on the AB stacked hexagonal lattice with the nearest- 
neighbor in-plane transfer integral 70 and the nearest- 
neighbor vertical coupling ji. 27 ~ 29 > These parameters are 
estimated as 70 « 2.8 eV and 71 w 0.4 eV. 19) 

Let iPai(Ra) (V'bz(-Rb)) be the amplitude of the wave 
function at R A (Rb) on the Ith layer (1 < I < TV), 
where Ra (Rb) represents the coordinate of an arbitrary 
site which belongs to the A (B) sublattice. Low energy 
states appear near the K + and K_ points in the two- 
dimensional Brillouin zone. The wave vector correspond- 
ing to the K± point is given by ±K = ±(2n/a)(2/3, 0). 
We express Vaz(-Ra) and V'bz(-Rb) as 



Vaz(Ha) = e iK - R -F+(R A ) + (-ly-V^i^A), 



MRb) = b F + (jRb ) + (_i)'e-'*-« B F m (R B ), 

(4) 

where F^(r) and F^(r) are envelop functions. Let us 
define the electron state \ i f?±) e near the K± point by 



N 



N 



l*±> e =E F A ; W|A ; ) e + E^W|Bz)e, (5) 
1=1 1=1 

and its time-reversed hole state by 

N 



\*±) h = j2(- 1 ) l - lF L(ry\A l ) l 



1=1 



N 



^(-iy^(rr|B ; ),. (6) 



Here |Aj) e (|B;) e ) represents the basis vector for electron 
states on the A (B) sublattice of layer I, and |A;)^ (|B;)/j) 
represents that for hole states. 

Within the effective mass approximation, one can 
show that these states satisfy H±\^±) e = e\^±) e and 
H±\&±)h = —e\ i &±)h, respectively, where H± is an ef- 
fective Hamiltonian given below 24 ) and e is the energy 
measured from \i. If the following basis set |Ai) p , |Bi) p , 
. . . , |Ajv) p , |Bjv) p (p — e or h) is chosen, H + is expressed 



as 



24) 
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V 



(7) 



J 



where Hi is the low-energy effective Hamiltonian for a 
monolayer graphene sheet 29 ) and V describes the inter- 



layer coupling. Here Hi and V are 

—p, 7^_ 



Hi 



V = 



7*4 






7i 




(8) 
(9) 



rj 

(3) ^ 



BdG 



A e ff(x)5 2 ATx27V 

-H+ 



(10) 



where 7 = (^/3/2)-f a and k± = k x ± ik y with k x — 
—\d x and k y = —id y . The matrix representation of H- 
is simply obtained by replacing k± — > k^ in H+. 

In the presence of superconducting proximity effect, 
one must treat the electron and hole states on the same 
footing, taking their coupling into account. For that we 
focus on the electron-hole space spanned by \®+) e and 
A natural, and certainly a possible way to proceed 
is to employ a Bogoliubov-de Gennes equation in the 
basis set adapted for the A^-layer SGS junction, |Ai) e , 
|Bi) e , . . . , |Ajv)ei |Bjv)e; \A-l)h, |Bi)fc, . . . , \A N ) h , \B N ) h . 
This has been indeed a conventional approach adopted 
for this system by several authors. 6 ) To be explicit, the 
Bogoliubov-de Gennes Hamiltonian in this conventional 
approach reads, 

A e ff {x)*52Nx2N 

where <5 2 jvx2JV = diag(l, 1,0, . . . , 0). Here A e g(x) is the 
effective pair potential, and S2NX2N represents the fact 
that only the first (top) layer is directly coupled with the 
superconductors. 30 - 1 This conventional approach, how- 
ever, has a drawback that it is not appropriate for study- 
ing the T- and T-depcndences of the Josephson effect 
taking a proper account of the planar structure of the 
junction. This point has already been mentioned in §1. 

Here, instead of naively applying the Bogoliubov-de 
Gennes equation we employ the tunneling Hamiltonian 
model proposed by McMillan. 25 ) This approach is espe- 
cially suited for a proper description of the coupling of 
the graphene sheet and the superconductors. The central 
theoretical tool of this approach is the thermal Green's 
function GJr,r';uj) with the Matsubara frequency lo — 
(2n+l)7rT. In this framework, the proximity effect is me- 
diated by quasiparticle tunneling between the graphene 
sheet and the superconductors, and is described by a 
self-energy. The thermal Green's function obeys 

{iujt z -H-Y,) G{r, r'; u) = T°S(r - r'), (11) 

where Z° _ = diag(l 2A rx 2A r, T 2WX 2iv) and t z = 
(l2JVx2JV,-l2JVx2Jv) with l 2 ATx2Ar being the 2N x 2N 
unit matrix, and H_ = diag(7T + , H + ). The self-energy E 
is represented as 20 ) 

£ _ / LoS 2 Nx2N A{x)&2Nx2N 

~ VA 2 + LO 2 \ A(x)*(5 2 7Vx27V -u52Nx2N 



x - 



(12) 



where F represents the strength of the tunnel coupling, 
and 6{x) is Hcaviside step function. The off-diagonal ele- 
ments are regarded as an energy-dependent effective pair 
potential, while the diagonal elements describe renormal- 
ization of a quasiparticle energy. If the w-dependence is 
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Fourier transformation: 



ignored in the self-energy by setting uj — 0, our model is 
reduced to the conventional one with the effective pair 
potential T. 6 ) 

3. Monolayer and Bilayer Cases 

We derive a general formula for the Josephson cur- 
rent in the monolayer and bilayer cases by using an 
analytical expression of the thermal Green's function. 
We employ the quasiclassical approach 21 ' 22 ) and focus 
on the slowly varying part of the Green's function, dis- 
carding the fast oscillating components on the order 
of the Fermi wave number. This allows for obtaining 
the Green's function analytically retaining sufficient ac- 
curacy at distances much longer than the Fermi wave 
length. 

3.1 Quasiclassical approximation 

By its construction, the Green's function G(r,r';co) 
inherits the matrix nature of the tight-binding (effective 
mass) Hamiltonian. In the case of monolayer graphene, it 
takes a 4 x 4 matrix form reflecting the sublattice and the 
electron-hole degrees of freedom. In the case of bilayer, 
the size of the matrix is further doubled (8 x 8) by the 
layer degree of freedom. The Green's function G(r, r'; to) 
in the bilayer case is defined in the electron-hole space 
spanned by |Ai) e , |Bi) e , |A 2 ) e , |B 2 ) e , \A-\)h, \^i)h, 
\&2)h, \^>2)h- However, within the low-energy regime of 
|e| <C 7i, we apply a second order perturbation theory 23 ) 
to the 8x8 Green's function and reduce it to a 4 x 4 func- 
tion defined in the reduced electron-hole space spanned 
by |Bi) e , |A 2 ) e , |Bi)/j, \A-2)h- This allows us to treat the 
mono- and bi-layer cases in parallel; the Green's function 
in the two cases are both represented by a 4 x 4 matrix 
form, which we denote as Gj{r,r';uj). The subscript j 
represent the number of layers: j = 1,2. The Green's 
function Gj(r 7 r';uj) obeys 

(iw7? x4 - Hj - £,) G,(r, p' ; W ) = f 4 ° x4( 5(r - r'), (13) 

where f| x4 = diag(l, 1, -1, -1), f 4 ° x4 = diag(l, 1, 1, 1), 
and Hj = dmg(H ■, H 3 ■) . The 2x2 effective Hamiltonian Although both the conduction and valence bands con- 
Hj is given in eq. (8) in the case of monolayer, whereas tribute to G jaa > , the latter is irrelevant in considering 



G j (x,x';q,u) = j d(y - j/) e -^v-v')Cf.(r, r'; u). (16) 

The fast spatial oscillations of Gj in the x direction 
are characterized by the wave number k, which is given 
by k = \J {p/"f) 2 — q 2 for the monolayer, and by k — 
x/fi/a — q 2 for the bilayer cases. Note that this value 
differs in the covered and uncovered regions due to the 
spatial dependence of ft. For later convenience we intro- 
duce the phase x<? = arg{fc + iq} in terms of which k 
and q are expressed as k = kp cos \q an d q ~ kp sin \q , 
respectively, where kp = /t/7 for the monolayer case and 
kp = \Jlija for the bilayer case. We separate fast oscil- 
lations of Gj by expressing it as 31 ) 

G j (x,x';q,cj)= £ G jaa . (x, x'; q, c)c ifc (— '*'). 

ct.ct'— ± 

(17) 

Within the accuracy of the quasiclassical approximation, 
the Green's function Gj a ^> obeys 

(iwf| x4 - Uj - tj) Gj aa * (x, x'; q, lu) 

= 5 a ^5(x-x')f° M , (18) 

where the 4x4 Hamiltonian Hj is given by Hj = 
dia,g(Hj + hj,Hj + hj) with 



hi 



—p, j(crk — iq) 

j(crk + iq) —ft 



jk x 
-fk x 



and 
H 2 



—a(ak — iq) 2 


—2a(ak — iq)k x 



-a(ak + iq) 2 



-2a(ak + iq)k x 



(19) 
(20) 

(21) 
(22) 



in the case of bilayer it takes the following form: 23 ) 



H 2 = 



-ak 2 _ 



-ak\ 



(14) 



where a = 7 2 /7i- The self-energy Sj is given by 20 ) 
f -ir ( co5^ 2 A(x)5% 2 \ 

2x2 



l — zxz — 

VA 2 + uj 2 \ A(x)*4x2 



low-energy properties since (i > A. Therefore we ex- 
clude the contribution from the valence band by using a 
unitary transformation which diagonalizes Hj . A pair of 
the eigenvectors, ipj a (q) and V'J (T (<?), of Hj are given by 

1 / e' iax " /2 



KM) 



x 9 



(15) 



1 

V2 



_ ae ia Xq /2 



(23) 



where 5. 



(i) 



(2~) 

diag(l, 1), and <5 2x2 = diag(l, 0). The ma- 



trix form of reflects the fact that only the top layer 
is in contact with superconductors in the bilayer case. 

Hereafter, we restrict our attention to the regime of 
moderate doping: 70,71 S> £t + U > \i » A. Assuming 
that our system is translationally invariant in the y di- 
rection under the condition of W 3> L, we perform the 



1 / e ^X q 

1 / Q^Xq 

KM = e-i'x. 



(24) 



and their eigenvalues are and — 2/i, respectively. Hence 
Hj is diagonalized as Uj a (qy HjUj a (q) = diag (0, — 2/i) 
with Uj a (q) = (ip^ a (q),'ipj r7 (q)), where the eigenvalue 
—2/i corresponds to the valence band. To exclude the 



4 



J. Phys. Soc. Jpn. 



DRAFT 



irrelevant contribution from the valence band, we per- 
form the transformation Gj ar7 > — Uj a (q)Gj (7cr 'iJj ( j'(q) 
with Uj a {q) = di&g(uj <7 (q),Uj a (q)) and retain only its 
(1,1)-, (1,3)-, (3,1)-, and (3, 3)-elements. 20 ' 30 ) Accord- 
ingly, we define G jar7 > by 



G 



jaa' 



[Gjacr'] 11 
3,1 



[Gj<7<j>] 13 
[Sjaa'] 3 3 



(25) 



which approximately satisfies the equation of motion for 
the quasiclassical Green's function 21,22 ' 

[iwr| x2 + iav jx d x T$ x2 - Y,j] G jaa > {x, x; q,w) 

= S a ^S(x - x')t 2 ° x2 , (26) 

where t| x2 = diag(l,-l), r 2 ° x2 = diag(l,l), and v jx is 
the x-component of the velocity given by Vj x — Vj-p cos Xq 
with t>iF = 7 and u 2 p = 2y/ajl. The self-energy is 



-iT,- 



VA 2 + w 2 



w A(x) 
A(x)* -w 



(9 



(27) 



where ri = T and T 2 = T/2. Note that E 2 is smaller by 
a factor of two than Ei reflecting the fact that only the 
top layer is in contact with the superconductors. 

Within the approximations described above, Gj aa > is 
expressed as Gj aa i — Gj aa > ® Aj aa >, where Aj arj i = 
Uj !7 (q)AUj r7 (q) with A = diag(l,0). The explicit forms 
of Aj a<7 i are 



Ai±± 
Ai 



1 



±e ±ix, 



1 

2 
1 

for the monolayer case, and 



A 2±± = \ 



A 



2±T 



1 

_ e Ti2x, 

e ±i2x, 
-1 



1 

Tl 



_ e ±i2x, 

1 



(28) 



-1 

3 Ti2x, 



(29) 



for the bilayer case. Consequently, Gj is represented as 
follows: 



Gj(x, x'; q, u) 

= ^2 G 3<r<r'( X > X '> < l> UJ ) 



Ajcrcr'C 



k(ax- 



(30) 



<T,cr'=± 



It is convenient to decompose Gj into four 2x2 Green's 
functions as 



Gj(x 7 x';q,Ld) 



gj(x,x';q,u>) fj(x,x';q,u>) 
f](x,x';q,uj) -gj(x,x';q,u>) 



(31) 

To obtain a formula for the Josephson current, we need 
only gj and fj. Using eqs. (26) and (30) we can show 
that in the uncovered region of — L/2 < x < L/2, they 
are expressed as follows: 

9j(x,x';q,u) 



— u~ 1 p ik ( x ~ x ') A 



+ v J^ k+ i*+*') Kj+ _ Cj+ _ 

+ v^e-^-^Xj— \-\e{x' -x) + Cj __] , (32) 



fUx,x';q,u>) 



= v ^ e ^^- k+ ^Aj ++ dj ++ 

+ vr x 1 e-^- x+k+ ^A j _ + dj. + 
+ v^e- i ( k -*- k+ ^A j __d j -, 



(33) 



where k^ — k±iu/vj X} and Cj aa i and dj aa i are unknown 
coefficients. 

3.2 General formula for the Josephson current 

We show that the Josephson current Ij ((f) is expressed 
in terms of Cj ++ and Cj__, and then determine these 
coefficients by applying a boundary condition at x = 
±L/2 to gj and /J. We finally derive a general formula 
for the Josephson current. 

The Josephson current is formally expressed as 

Ij(<p)=4W I ^-Tj2tv{j jx9j (x;q,uj)}, (34) 

w 

where the factor 4 comes from the spin and valley degen- 
eracies, the current operator Jj x is defined by 



Jix = ej 



1 

1 



•hx = -2ea 



k x +\q 



(35) 



(36) 



k x -iq 

and gj (x; q, u) = [g 3 (x, x - 0; q, w) + g 3 {x, x + 0; w)]/2. 
Substituting eq. (32) into eq. (34) and changing the in- 
tegration variable from q to x 9 , we obtain 

I j (<p)=2eN j (0)W 
/■+f 

x / ^ dx g Vj x Tj2(c j+ +(v)-Cj—(<p)), (37) 

where Nj(0) the density of states per spin including 
the valley degeneracy is given by Ni(0) = fi/(ir-f 2 ) and 
N 2 (0) = 1/(2™). 

The unknown coefficients arc determined by the 
boundary condition at x = ±L/2 for gj(x,x';q,uj) and 
fj(x,x';q,uj), 32 ) which we outline below. It is necessary 
to distinguish Xq m the covered and uncovered regions 
in the following argument, so we rewrite Xq as (<f>) i n 
the covered (uncovered) region: 



Xq 



{-L/2 <x< L/2) 
(L/2 < \x\) 



(38) 



To represent the boundary condition we introduce elec- 
tron wave function \I>j^ (x; q, oj) and hole wave function 
ty^j(x;q,uj) in the uncovered region. They are decom- 



-it/ (x 



u 3++\ 
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posed into right-going and left-going components as 

^j(x;q,u) = c N+ (x,u)ip c j+ (q) + c N _(x, w)^ c _(g), 



which are 



(39) 



^j(x;q,u) = d N+ (x,u})tf + (q) + d N -(x,u})tf_(q), 



where c^±(x, co) (d^±(x,co)) contains the x-dependent 
factor of e ±lk+x (e ±lk x ). Applying the prescription by 
Zaitsev, 31 ) we derive a relation between cn±(±£/2,w) 
and d^±(±L/2,u>). We leave its derivation to Appendix 
A, and here refers only to the final result, 

cn+(±L/2,w) \ / d N+ (±L/2,u>) ' 

c N _(±L/2,u;) J-^^/^l d N -(±L/2,u) 



(41) 



with 



Bj(±L/2) 



where 



5 ±i V /2 

A 



r 3 



A = 



Vw 2 + A 2 



A, 



V^ 2 + A 2 ' 



o = Vc^ + A 2 , 



(42) 

(43) 
(44) 
(45) 



and Tj and IZj = 1 — 7} denote the transmission and 
reflection probabilities for an electron at the Fermi level 
being incident at the uncovered-to-covered interface with 
an incident angle </>, and £ is the phase of the correspond- 
ing reflection coefficient. The transmission probability for 
the monolayer case is 

^_ 2 cos <j> cos 9 



1 + cos(0 + 6»)' 
and that for the bilayer case is 

4 (1 + C) 2 w cos <j> cos 9 



r 2 = 



(l + £ 2 )A + 2£B 



(46) 



(47) 



with 



(48) 
(49) 



A = 1 + w 2 + 2w cos(0 + 9), 
B = cos 2(p + w 2 cos 29 + 2w cos(0 - 9) , 
where w = \J Ji/Jp + U) and 

{^J 1 + sin 2 <j) + sin 0) ( + sin 2 <j> + U/n — sin< 

£ = — — — 7 

+ sin 2 — sin </>^ ( ^/l + sin 2 <f> + U / + sin < 

(50) 

The derivation of eq. (46) and (47) is given in Appendix 
B. We determine the unknown coefficients in eqs. (32) 
and (33) by using eq. (41) as the boundary condition. 
Equation (41) provides us coupled eight relations, two of 



(-i + c J++ )c : 
c,_ + c- ife+i / 2 



fc+L/2 



= Bj(L/2) 



a J++ c 



ifc^L/2 I 

(51) 



(40) Solving these coupled relations we obtain 



with 



W^) = 4(J- 1 ) (^) 



Tj = (lu 2 + _4 2 £! 2 ) sinh Kj + 2AjU)Cl cosh Kj 

+ iB 2 n 2 sm2kL + iA 2 smtp, (53) 

IIj = (Cj 2 + A 2 Q 2 ) cosh Kj + 2AjU)£l sinh Kj 

-B 2 n 2 cos2kL + A 2 cos 93, (54) 

where Aj = (1 + TZj)/Tj, Bj = 2^/TZj/Tj, and 
Kj = 2u>L/vj x . The other coefficient Cj—{ip) satisfies 
Cj--(<p) = Cj ++ (—ip), so we do not present it explicitly. 

Substituting these results into eq. (37), we finally ob- 
tain a general formula for the Josephson current 



I j { i p) = eAf j 



> cos 4>T 



' sin ip 



IT 



(55) 



where Afj = 2vj-pNj(0)W represents the number of con- 
ducting channels. Using this general formula one can nu- 
merically calculate the Josephson current in the planar 
junction for arbitrary parameters. If we take the strong- 
coupling limit of Tj 3> Ao (Ao: the magnitude of the pair 
potential at T = 0), eq. (55) is reduced to the expression 
for the Josephson current in an SNS junction (N: nor- 
mal metal) with same barriers at the two NS interfaces, 
derived by Galaktionov and Zaikin. 32 ) 

4. Extension to the TV-layer Case 

To study the Josephson effect in a planar junction of 
A-layer graphene, we first decompose the 4N x 47V ther- 
mal Green's function into a set of 4 x 4 Green's functions. 
With this decomposition we can easily derive a formula 
for the Josephson current by applying the argument pre- 
sented in the previous section. 

Following Koshino and Ando, 24 ) we decompose H + for 
an A-layer graphene sheet into 2x2 monolayer-type 
and 4x4 bilayer-type Hamiltonians with the help of 
an appropriate choice of the bases. For N being odd, 
the 2N x 2N Hamiltonian H + is decomposed into one 
monolayer-type Hamiltonian and (N — l)/2 bilayer-type 
Hamiltonians. While, for N being even, H + is decom- 
posed into A/2 bilayer-type Hamiltonians. This proce- 
dure can be adapted to a multilayer graphene sheet in 
planar contact with superconductors as demonstrated in 
ref. 30. 

Applying this decomposition also to the thermal 
Green's function G(r,r';uj), we can decompose it into 
4x4 monolayer-type and 8x8 bilayer-type Green's func- 
tions. The obtained 8x8 bilayer-type Green's function 
is then reduced to a 4 x 4 Green's function. These de- 
composition and reduction procedures are explained in 
Appendix C. In the case of N being odd, G(r,r';co) is 
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decomposed into one monolayer-type Green's function 
GN,o(r,r';uj) and (N — l)/2 bilayer-type Green's func- 
tions GN iTn (r, r'; u>) with m = 2, 4, . . . , N — 1. While, 
for TV being even, G(r,r';w) is decomposed into 7V/2 
bilayer-type Green's functions G n , r' ] u>) with m = 
1,3,...,/V-1. 

The Green's function GM,m obeys 

(iwf| x4 - Fjv.m ~ Sjv,m) <5jv, m (r,r';<*>) = f" x4 c5(r - r'). 

(56) 

For the Green's function of monolayer-type with m = 0, 
the Hamiltonian and the self-energy are given by Hn,o — 
diag(i?i, H\) and 



-ir 



AT,0 



A(^)4x } 2 



VA 2 + u; 2 \ A(a;)*4 x2 -wA 



'2x2 



x 6 ( |a;| - - 



with 



Tat ,o 



(57) 



(58) 



7V + 1 

For the bilayer-type with m/0, the Hamiltonian is given 

jN.m ttJV 
»2 '""2 



by Fjv >m = diag( J ff 2 A ' ,m , J ff^ m ), where 



if. 



JV,m 



^/fc 2 



-A 



with 



A«, m = 2 sin 



77J7T 



2(JV + 1)7 ' 



(59) 



(60) 



The self-energy is given by 

_ — 2irjv !m / 
N ' m ~ VA 2 + oj 2 \ A(x)*5. 



(2) 
2x2 

* X ( 2 ) 
2x2 



-ljS, 



(2) 
2x2 



with 



x I M - - 



2 cos 2 



N,m 



2(JV+1) 



7V + 1 



r. 



(62) 



The Green's function Gjv,m is reduced to that in the 
monolayer (bilayer) case if N = 1 (TV = 2). 

Applying the procedure described in the previous sec- 
tion to Gjv.m; we obtain the corresponding contribution 
lN,m(<P) to the Joscphson current, which is given by 
cq. (55) with the following replacements: 



N,m 

VjF ->■ V F , 

Nj(0) ^ N N , m (0), 

Fj — > rjv,m: 



where 



JV.O 

% ' =Vif = 7, 



N.m 



V2F 



an 



(63) 
(64) 
(65) 

(66) 
(67) 



and 



JVjv,o(0) = JVi(0) = -^, 

iVjV,m(0) = \ N ,mN 2 (0) 



A 



iV., 



27ra 



(68) 



(69) 



The total Josephson current In{<p) m the TV- layer case 
is expressed as 

M¥>) = £jiV,m(¥>)- ( 7 °) 
m 

5. Analytical and Numerical Results 

Let us focus on the short-junction limit: L <C £, where 
£ = WjF/(27rAo) is the superconducting coherence length. 
This limit has a particular importance to the interpreta- 
tion of experiments. Taking this limit helps simplifying 
the denominator of eq. (55) by setting cosh/tj w 1 and 
sinh Kj w 0, resulting in 

f+f 

(09) = eA/j / d0 cos ^ 
A 2 sin 09 



r, 



ft 2 



75 



fi 2 cos 2fcL + A 2 cos ip 
(71) 

Based on this formula, we first study analytically the 
behavior of the Josephson current in several limiting 
cases, then estimate the expression numerically to ob- 
serve those features that are not captured by the analyt- 
ical arguments. 

5.1 Analytical result 

5.1.1 Cases of homogeneous carrier density 

Let us first consider the simplest case of U = (with- 
out carrier inhomogeneity) . Since Tj = 1 and IZj = in 



(61) this case, eq. (71) is further simplified to 



I j (<p)=2eAf j Tj2- 



A sin w 



Cj 2 + n 2 + A 2 cos 09 



(72) 



for the monolayer and bilayer cases. 20 ) The behavior 
of Ij(ip) predicted from cq. (72) has been discussed 
in ref. 20. In the strong coupling limit of Tj 3> Ao, 
this formula reproduces the result derived by Kulik and 
Omel'yanchek (KO), 33 ) 



09 /A cos £ 
Ij (09) = eA/j A sin - tanh I ^ 



(73) 



In the weak coupling limit of Ao 3> Tj, the Josephson 
current behaves as 



Vx—i. I r i cos 2 



Ij(ip) = eMj Tj sin - tanh > ^ 



in the low-temperature regime of Tj T, 34 ) and 

r 2 A 2 

Ij(<p) = eA/j sin 09 



(74) 



(75) 



in the high temperature regime of T c > T. Note that 
the T-dcpcndence of Ij (09) in the weak coupling limit 
is qualitatively differ from the KO result which yields 
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I(ip) = ej\fA 2 {AT)~ x sin 95 near T c . The behavior re- 
vealed in cq. (75) has not been reported in literature, 
indicating that the Josephson effect in a planar junction 
is significantly affected by the coupling strength. 

5.1.2 Back to the generics — cases of inhomogeneous 
carrier density 
Let us go back to eq. (71), and investigate the behavior 
of Ij(<p) in several specific limits. Rewriting A 2 /IX,- in 
cq. (71) as 



(r,A) 2 



II, $p j -2(r j A) 2 sin 2 f 



with 

% = 1 + 



r 3 



cos2fc7, 



(76) 



(77) 



Vj =cu 4 + (a 2 + 2r jV / A 2 + ^ 2 + r 2 ) lo 2 + (r,A) 2 , 

(78) 



we approximate 2Tj\/ A 2 + lo 2 in T>j by 2TjA. This is 
justified in three situations: the low-temperature limit 
of T c 3> T, the strong-coupling limit of Tj ^> A , and 
the weak-coupling limit of Ao 3> Tj. Under this approx- 
imation we perform the summation over w in eq. (71) in 
terms of a contour integral, and obtain 



Ij(<p) = ej\fj(TjA) 2 simp 
tanh 



d(j) 



COS( 



Ej-L —Ej 



4T 



tanh 



E j + +E 3 



4T 



E j+ ~ E 3- 



E 3 + + E 3 



where 



with 



^• ± = A /(A + r j ) 2 ±2 v '|/ j |r j A 



f _ 1 _ 2sin2 f 



(79) 



(80) 



(81) 



Equation (79) reduces to eq. (31) of ref. 20 in the homo- 
geneous case of U = O. 35 ) 

The characteristic behavior of Ij ((p) in the strong- and 
weak-coupling limits follows from eq. (79). The Joseph- 
son current in the strong-coupling limit is given by 



Ij(ip) = ejVj A sin <p 



dcj) 



COS ( 



■. tanh 



2T 



(82) 



for an arbitrary T. In the weak-coupling limit, we obtain 

cos<£ . , / \/\fi\ T A 



r+i 
r jTj sirup j 



d(j) 



: tanh 



2T 



in the low-temperature regime of Tj 3> T, and 



Cj(TjA) 2 simp 
48T3 



(83) 



(84) 



in the high-temperature regime of T c > T, where 



d(j)- 



cos< 



(85) 



Note that eqs. (82), (83), and (84) reduce to cqs. (73)- 
(75), respectively, in the homogeneous case of U = 0. 

5.1.3 Zero-temperature behaviors 

Let us consider the critical current Ij c at T = 0, where 
cq. (79) simplifies to 



Ij(<p) = ejVj 



A + r, 



sin ip 



cos< 



2$, 



1 + 



2r,A 



In the homogeneous case of U = 0, <£,- 



2 since Tj 



(86) 



1 irrespective of (j). In this case Ij(ip) is maximized at 
<p = ir (mod 2n). Thus the critical current at U = is 
determined as lf=° = eMjTj A /(A + Tj). This yields 

2fiW TAq 



T u=o 

1 lc 



A + r 



for the monolayer, and 



T u=o 

1 2c 



2fiW 

7T7 1 





§A 




Ao + i 



(87) 



(88) 



for the bilayer cases. Speaking of the order of magnitude, 
one can notify that the critical current in the case of 
bilayer is greater than that of monolayer by a factor of 

5.1.4 Cases of multi-layers 

The critical current for the cases of multi-layers at T = 
can be discussed on the same footing. As demonstrated 
in eq. (70), the Josephson current in the case of iV-layer 

(1) is given by the summation of ijv, m (v?) over m, 

(2) and every contribution is maximized at ip = ir when 
17 = 0. ' 

Combining these two observations, we see that the crit- 
ical current for the iV-layer case is simply given 
b y E m lN,m(<P)\ v ^, where I N , m {p>) is obtained by per- 
forming the replacements, eqs. (63)-(65), in cq. (86). 
Listing up the cases of small N, one can find 



(89) 




T u=o 



(90) 



for the tetra-layer cases. 
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7=2 (bilayer case) 



10 



U/\L 



20 



Fig. 2. The normalized critical current Ijc/lf c ^° at T = as a 
function of U/fi, where r = 0.1 (dashed lines), 1.0 (dotted lines), 
10.0 (solid lines). The data for r = 0.1 and r = 10.0 arc almost 
overlapped and indistinguishable in both the monolayer and bilayer 




77 r„ 



Fig. 3. The normalized critical current Ijc/ljf ^ at U/fj, = as 
a function of T/T c . This result applies to both the monolayer and 
bilayer cases. 



5.2 Numerical result 

The use of the analytic formulas obtained so far is 
limited in the vicinity of a particular limit we have spec- 
ified each time in the space of control parameters, (U, T). 
Here, to give an insight on the behavior of the criti- 
cal current in the entire range of this parameter space, 
we discuss some numerical plots obtained by directly 
estimating the expressions, such as eqs. (86), (72) and 
(71) (instead of trying to simplify them by restricting 
the range of validity). In the actual computation, the 
following set of parameters is employed: 70 = 2.8 eV, 
71 = 0.4 eV, a = 0.246 nm, A = 120 ,ucV, \i = 32 mcV, 
and L = 200 nm. 

5.2.1 V '-dependence of the critical current 

Fixing the temperature at T = 0, let us focus on 
the [/-dependence of the critical current in the mono- 
layer and bilayer cases. At T = the critical current Ij c 
can be computed at an arbitrary value of U, measure 
of the carrier inhomogeneity, by numerically estimating 
eq. (86). In Fig. 2 the critical current Ij c (normalized by 
Ijlr ) is shown as a function of U//J, for different val- 
ues of a parameter r = Tj/Ao that characterizes the 
coupling strength between the graphene sheet and the 
superconductors. In the figure, different curves, corre- 
sponding respectively to a different choice of this param- 
eter r = 0.1, 1.0, and 10.0, and to a different number of 
layers (j = 1, 2) are superposed for comparison. The be- 
havior of Ijc/Ijir shown in the figure indicates that it is 
almost insensitive to the variation of r in this range, but 
exhibits a qualitatively different behavior in the mono- 
layer and bilayer cases. 

The suppression of ij c //j£ =0 as the increase of U is 
more pronounced in the bilayer case. In the case of mono- 
layer, the ratio ii c /-/j^ =0 , indeed converges to a constant 
value in the limit of U — > 00. 6 ' On contrary, hc/I^ 
decreases monotonically toward zero. This contrasting 
behavior reflects the difference of the transmission proba- 
bility Tj in the two cases. In the simple case of an electron 



at the Fermi level incident perpendicularly to the inter- 
face, the electron is completely transmitted in monolayer 
even though U is very large. Contrastingly, the transmis- 
sion probability significantly decreases with increasing U 
in the bilayer case. 

5.2.2 T-dependence of the critical current 

The T-dependence of the critical current is fully en- 
coded in eq. (71). Here, we estimate this formula nu- 
merically, focusing on the cases of the monolayer and bi- 
layer systems. The actual computation is performed for 
r = 0.2, 1.0, 5.0, and r — » 00. The amplitude of the pair 
potential is determined by the gap equation 



l = Ai, 



L 



de tanh ( — 6 

\ 2T 



/v^ + A 2 , (91) 



where A; n t is the dimensionless interaction constant, and 
the Debye energy is chosen as ed/Aq = 200. For com- 
parison we have considered the cases of homogeneous 
(U/fi — 0) and highly inhomogeneous (U//J, = 10) car- 
rier density. 

In Fig. 3, the critical current Ij c (normalized by 1^ = 
&A/}Ao) is shown as a function of T/T c for the homo- 
geneous case of U/fi — 0. Recall that at U/fi = 0, the 
Josephson current for the mono- and bi-layer systems is 
described by the single formula of eq. (72). Thereby, the 
T-depcndcnt normalized critical current for j = 1 and 
j = 2 obviously coincides in this limit. 

We observe that Ij c is a concave function of T in the 
strong-coupling limit of r — > oo, where the KO result is 
reproduced. However, it crosses over to a convex function 
with decreasing r. Such a convex T-dependence was not 
observed in the previous study 11 -* based on an energy- 
independent effective pair potential model. The coupling 
strength Tj crucially affects the T-dependence of the crit- 
ical current. 

The behavior of the critical current in the highly in- 
homogeneous case of U I [i = 10 is shown in Fig. 4. Apart 
from the saturation of the normalized critical current 
at low temperatures (and this happens irrespective of 
the value of r), it seems that an overall character of 
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7=2 (bilayer case) 
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Fig. 4. The normalized critical current Ij c /I^ at Uf^i = 10 for 
the (a) monolayer and (b) bilayer cases as a function of T/T c . 



the T-dcpendence of the critical current remains un- 
changed both in the monolayer and bilayer cases. This 
leads us to the following statement: "the concave-to- 
convex crossover of the critical current, which occurs as 
the decrease of the coupling strength, which also appears 
irrespective of the homogeneity of the carrier density, 
should be regarded as a characteristic feature of the pla- 
nar Josephson junction" . 

6. Summary 

After a detailed formulation of our model proposed 
for the planar Josephson junction of graphene, we have 
derived and estimated a class of formulas for the Joseph- 
son current. The derivation was done in the framework 
of quasiclassical approximation, and the formulation was 
generalized to be applicable to the case of multilayer 
graphene. The obtained formulas allow us to estimate the 
Josephson current numerically in the planar junction of 
A^-layer graphene, for an arbitrary choice of the set of pa- 
rameters, that includes T (temperature), L (separation 
of two superconductors), V (coupling strength between 
a graphene sheet and superconductors) , and U (effective 
potential controlling carrier inhomogeneity) . Throughout 
the paper, we have assumed that the chemical potential 
\i is away from the Dirac point. 

Much emphasis has been put on the behavior of the 
Josephson current in the short-junction limit in the 
monolayer and bilayer cases. It was demonstrated that 
the coupling strength crucially affects the temperature 
dependence of the critical current in an unexpected man- 



ner. This should be regarded as a characteristic feature of 
the planar junction. We have also shown that the depen- 
dences of the critical current on \i and U differ qualita- 
tively in the monolayer and bilayer cases. The difference 
in the [/-dependence (see Fig. 2) reflects the contrast- 
ing behavior of the transmission probability across the 
uncovered-to-covered interface in the two cases. 
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Appendix A: Derivation of the Boundary Con- 
dition 

The unknown coefficients Cj ++ (ip) and Cj (ip) in the 

general formula for the Josephson current [eq. (37)] have 
been fixed by the boundary condition [eq. (41)] at x = 
±L/2. Here, focusing on the case of x = L/2, we give 
an explicit derivation of eq. (41). The derivation consists 
of three steps. Let us start with electron and hole wave 
functions in the uncovered region 

^j{x;q,u}) = c N+ (x,u)rf + (q) + c N _ (x, w)^ c _ (g), 

(Al) 

*Nj {x;q,u)= d N+ (x, u)i>j + (q) +d N _(x, uj)^_ (q) , 

(A2) 

and those in the covered region 

^iy fa; = cs+(x,u})ip] + (q) +c s _(x,w)^ c _(<?), 

(A3) 

*|j (x; q, u) = d s+ (x, uj)^ + (q) + d s _ (x, w)^ c _ (q), 

(A4) 

where cn±(x,u>) and cs±(x,cj) contain the x-dependent 
factor of e x and d-^±(x, ui) and ds±(x,uj) contain 
the factor of c ±lfe x . Note that il>j±(q) and k depend on 
whether x is in the covered or uncovered region. Firstly 
we relates cs±(L/2,lj) and ds±(L/2,±) on the basis of 
the Bogoliubov-de Gennes equation. Secondly we derive 
the relation between cn±(£/2,w) and cs±(L/2,u>) and 
that between d-^±(L/2,w) and ds±(L/2,u>) following the 
prescription presented by Zaitsev. 31 ^ Finally we derive 
eq. (41) by combining these relations. 

In order to relate cs±(L/2,u>) and ds±(L/2,u>), we in- 
troduce the Bogoliubov-de Gennes equation 

[ iwT 2x2 + ivVj x d x T$ x2 - £,] ^sja(x; q,u>) = 0, (A-5) 

where a = + (— ) specifies the right-going (left-going) 
component. In the region of x > L/2, eq. (A-5) is simpli- 
fied to 

iuj + \avj x d x iAe 1<p / 2 \ T , . 

iAc-W2 i, ■ J «,«)=o. 

(A6) 



Its solution is 
$>Sjcr(x;q,uj) 



A , 

e-W2 



-(x-L/2) 



(A-7) 
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The corresponding four-component solution ^sja is ob- 
tained by adding the pseudospin (sublattice) degree of 
freedom in accordance with ipj a (q), leading to 



with 



^Sja{x;q,uj) 



-(x-L/2) 



(A- 



Comparing this with eqs. (A-3) and (A-4), we easily see 
that the electron and hole components satisfy 



c s+ (L/2,uj) 
cs-(L/2,u) 



A 







A 



d s+ (L/2,u) 



(A-9) 



We next derive the relation which connects 
c^ a {L/2,ui) and cs<r(L/2, w), and d^ a {L/2,uj) and 
dso-(L/2,iu). Since the effective potential drops to — U 
across the interface at x ~ L/2 within the length scale 
comparable to, or shorter than the Fermi wave length, 
we are allowed to describe wave functions in this region 
by ignoring the influence of lo and A. 31 ) Hence, with 
this approximation, we can connect c-^ a (L/2,uj) and 
cs<y(L/2,uj) by considering a scattering problem at the 
Fermi level ignoring the coupling between electron and 
hole components. Let tj and rj be the transmission 
and reflection coefficients for an electron incident at 
the interface from the x < L/2-side with an incident 
angle <p. In terms of tj and rj we can construct two 
independent wave functions: 

■ [c N+ {x)ip c j+ {q) 
ipij(x;q)={ +r 3 c N _(x)ip c J _(q) 
--t jCs+ (x)ij c 3+ (q) 



ip2j(x;q) 



/ cos 

T-^CN-O^-fa) 

7 J==[ Cs _(x)i>U<l) 



Or < L/2) , 
(x > L/2) 

(A-10) 



-r jCs+ (x)^ + (q)] 



(x < L/2) 

(x > L/2) 

(A-ll) 



Here t,- 



and fj are the transmission and reflection co- 



efficients for an electron incident from the right, and 
they are given by tj — tj and fj = ^ r jtj/tj- A gen- 
eral wave function is represented by their superposition 
as tjjj — aipij + bip2j ■ Comparing this with eqs. (A-l) and 
(A-3), we obtain 



c N +(L/2,w) \ = 1_ 

c N -(L/2,cj) J ^/cosl 

c s+ (L/2,lu) \ = I 

c s _(L/2,w) J \ 



These two equations yield 



I 



(A-12) 



(A-13) 



c N+ (L/2,o;) 
c N _(L/2,w) 



lcos9 M ( c S+ (L/2,w) 
COS0 j { c s _(L/2, W ) 



(A-14) 



Mj = 



j 

holds 



(**) 



3' 

for 



(A-15) 

d-^±(L/2,u>) and 



The same relation 
ds±(L/2,w). 

Eliminating cs±(L/2, co) and rfs±(L/2, lo) in cq. (A-14) 
and the corresponding relation for d^±(L/2,w) and 
ds±(L/2,cj) with use of eq. (A-9), we finally arrive at 
eq.(41). 

Appendix B: Transmission Probability across 
the Interface 

Here, we derive the expressions of the transmission 
probability for an electron at the Fermi level being in- 
cident at the interface at x = L/2 from the x < L/2-side 
with an incident angle 4>. For simplicity we shift the ori- 
gin of the one-dimensional coordinate to x = L/2 and 
replace x — L/2 — > x. 

Let us first consider the case of monolayer. Using the 
eigenfunctions of Hi , one can express the two-component 
wave function for this scattering problem as 



ipi{x;q) 



sjl COS 


V e^ 2 


+n 


( 


Sjl COS \ 


t eil "" 


/ e -ie/2 


1 V2 cos 


9 ( e i0/2 



e vl \ 
_e-W2 J ^ <0 ) ' 

(x > 0) 

(B-I) 

where the elements of ip1± (q) is explicitly shown for clar- 
ity. Note that the value of k = \J (A/7) 2 — Q 2 differs in 
the left and right of the interface. By matching the two 
components of the wave function at x = we obtain 



h = 



n 



^/cos< 



cos 



6 



sin ■ 



i±e 
2 

—6 



COS 



4>+0 



(B-2) 



(B-3) 



We easily see that the transmission probability 71 = \ti\ 2 
is given by eq. (46). 

We next consider the bilayer case in which evanescent 
modes play a role. 36 ' 37 - 1 By using the eigenfunctions of 
H 2 , the two-component wave function is expressed as 

-e- 1 " 



ip2(x;q) 



V2 cos 

+r 2 



^2 cos ^ 



e 

-e i( 



-ac" 



*2 



V2 cos 

+bc- 



1 

K — q 

n -i6 



K+q 



{x < 0) 



{x > 0) 

(B-4) 
^ Note 



where fc = y/ (ft/ a) — q 2 , and k = yj (ft/ a) + q* 
that the value of k and k differ in the left and right of 
the interface. It should be noted that the wave function 
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contains the probability amplitudes at Bi and A2 sites. 
Remember that the effective Hamiltonian H2 for the bi- 
layer case is obtained in the framework of a second order 
perturbation theory. Within this perturbation theory the 
probability amplitudes at Ai and B2 sites are obtained 
from those at Bi and A 2 sites. Indeed, for a given two- 
component wave function ip2(%] q) = t (u(x), v(x)) for J$i 
and A 2 sites, we can show that the corresponding four- 
component function is given by 38 ' 
t 

tp 2 (x;q) = 



7 ~ 

— k+v(x), u(x), v(x), 
7i 



- — fc_u(x) I , 

71 



(B-5) 



where k± = k x ±iq. Constructing tp2{x;q) from eq. (B-4) 
and then matching its each component at x = 0, we 
finally obtain 



*2 = 



r 2 = 



2 (1 + C) y/wy/cOS (j) cos 9 



c-^£(l + weH^+e)) + c^(l + w)c-^+ 6 ) ' 



(B-6) 



3^£(1 



we 



-K4>- 



»)) _ e -i*(i _ we^-V) 



e-^£(l + we i W+ )) + e^(l + we- W)) ' 



where w = \J + £/) and 



<? 2 + 9) 



£ = 



) h/^ 



q 2 - q 



(B-7) 



(B-8) 



Equations (47) and (50) are easily obtained from 
eqs. (B-6) and (B-8), respectively. 

Appendix C: AT- layer Thermal Green's Func- 
tion G N ,m 

Applying the method presented in ref. 30, we decom- 
pose the 47V x 47V thermal Green's function G(r,r';uj) 
into a set of 4 x 4 functions. Let us define 



TV + 1 
2 

VnTT 

2 



(-iy 



cos 



2(TV + 1) 



(C-l) 
(2i-l)), 



sin 



mir 



(C-2) 
-pi)), (C-3) 



y/N + T ' \2(N + 1)' 

for a given TV. In terms of these functions we define the 
following bases for the odd- TV case: 

|1 ) =<A)(l)|Al) e + ^o(2)|A 3 )e 

+ ■■■ + M(N + l)/2)\A N ) e , 

|2 > =Po(l)|Bi>e+W>(2)|B 3 >e 

+ --- + ^ ((TV + l)/2)|B w ) e , 
|3 ) =^o(l)|Ai) fc + ¥Jo(2)|A3) fc 

+ --- + <p ((N + l)/2)\A N ) h , 
|4 ) =W)(l)|Bi) fc + V o(2)|B 3 )fc 

+ --- + ¥>o((JV + l)/2)|B JV ) fc) 



(C-4) 
(C-5) 
(C-6) 
(C-7) 



and 

= V +(l)|Ai) e + V +(2)|A3)e 

+ --- + <p+((N + l)/2)\A N ) e , (C-8) 

= V +(l)|B 1 ) e + V +(2)|B 3 )e 

+ ... + ^+((TV + l)/2)|B JV > e , (C-9) 
= ^-(l)|A 2 ) e + ^-(2)|A 4 ) e 

+ --- + ^-((TV-l)/2)|A JV -i> e , (C-10) 

= V -(l)|B 2 ) e + V -(2)|B 4 )e 

+ --- + ^-((7V-l)/2)|B JV _ 1 ) e , (C-ll) 

-V+(mi)h + V+(2)\A 3 ) h 

+ ... + tp +((N + l)/2)\A N ) h , (C-12) 

= V +(l)|B 1 ) fc + V +(2)|B 3 ) fc 

+ --- + <p+((N + l)/2)\B N ) h , (C-13) 

= V -(l)|A 2 ) h + ¥)-(2)|A4)fc 

+ --- + < / >-((TV-l)/2)|A JV _ 1 ) /l , (C-14) 

= V -(l)|B 2 ) fc + V -(2)|B 4 ) fc 

+ --- + ^-((/V-l)/2)|B Ar _ 1 ) /l . (C-15) 

In the odd-TV case, we can decompose the 4TV x 4TV 
Hamiltonian H_ into one monolayer-type Hamiltonian 
with dimensions 4x4 and (TV — l)/2 bilayer-type Hamil- 
tonians with dimensions 8 x 8 by transforming the orig- 
inal bases to the new ones presented just above. 30 ) Ap- 
plying this procedure, the 4TV x 4TV thermal Green's 
function G(r 7 r';uj) is decomposed into one monolayer- 
type Green's function Gn,o{t, r'\ uj) and (TV — l)/2 
bilayer-type Green's functions GAr, m (r, r'; lo) with m = 
2,4,...,TV — 1. In doing so we approximately ignore weak 
coupling between different Green's functions due to the 
self-energy. 30 ) The 4x4 monolayer-type Green's function 
Gn,o is defined so that its (i, j)-element — 1,2,3,4) 
is given by 



[Gjv,o] , 



(io|G|j > 



(C-16) 



From eq. (11) we can show that this obeys eq. (56). The 
8x8 bilayer-type Green's function Gat,™ is defined so 
that its (i, j)-element (i, j = 1, 2, 3, . . . , 8) is given by 



G 



h3 



= (irn\G\j m ). 



(C-17) 



From eq. (11) we can show that this obeys 

H N , m - Sjv,m) G~N, m (r, r'\ uj) = xS 8(r - r') 

(C-18) 



3x8 



with 



8x8 



'2x2- 



^ixii ^8x8 

. The Hamiltonian H. 



r, 



N,m 



IS 



H. 



N.m 



Hi 



#1 



(C-19) 
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with \N, m given in eq. (60). The self-energy is 



-2ir 



'2x2 



CO 



2 



x2 



02x2 
02x2 



A(x) 

-co 







L 



x\ — 



(C-20) 



where T N>m = (l/2)<^+ (1) 2 T. Now we reduce the 8 x 8 
Green's function G^m to a 4 x 4 Green's function in 
accordance with the argument presented by MacCann 
and Fal'ko. 23 ) The interlayer coupling (i.e., 71 in V 
and V^) forms dimers in the electron space spanned by 
|l m ) and |4 m ), and in the hole space spanned by |5 m ) 
and |8 m ). Since the energy of dimer states is greater 
than 71, we can safely ignore these states as long as 
/i -C 71. Consequently, within a second order pertur- 
bation theory, low-energy quasiparticles are described in 
the electron-hole space spanned by |2 m ), |3 m ), |6 m ), |7 m ) 
and the corresponding 4x4 Hamiltonian is given by 
HN,m = diag(i?^ v,m , H2' 171 ) with the reduced bilayer- 
type Hamiltonian H 2 ' m presented in eq. (59). In accor- 
dance with this reduction, GN,m is reduced to the 4x4 
Green's function Gn,™, which obeys eq. (56). 

In the even-iV case, the 4N x 4N thermal Green's 
function is decomposed into N/2 bilayer-type Green's 
functions GN,m{r, r'; lo) with m = 1, 3, . . . , N — 1. The 
bilayer-type Green's function GN,m{r, r'; u) is defined so 
that its (i, j)-element (i, j = 1, 2, 3, . . . , 8) is given by 



G 



N,m 



(i m \G\j m ) 



with 



U =¥>+(!) 



2m) =¥>+(!) 



3m) = Pmi 1 ) 



4m) =<P m (l) 



5m) =vtl(X) 



6m) =^(1) 



7 m ) =<p-(l) 



8m) = ¥>"(!) 



Ai 
+ 
Bi 



B 2 
+ 
Ai 
+ 
Bi 
+ 
A 2 
+ 
B 2 
+ 



e + ^+(2)|A 3 ) e 

•• + ^+(7V/2)|A A r_ 1 ) e , 

e + ^+(2)|B 3 ) e 

■■ + Cf+(N/2)\B N - 1 ) e , 

e+^m(2)|A 4 ) e 

•• + ^ m (7V/2)|A A r) e , 

e +^ m (2)|B 4 )e 

■■ + V -(JV/2)|B JV ) e , 

fc + ^(2)|A3) fc 

■■ + V +(JV/2)|A JV _ 1 ) fc , 

fc+¥>+(2)|B 3 >/, 

■■ + V +(JV/2)|B JV _i) h) 

^ + ^m(2)|A 4 ) ?i 

■■ + <p-(N/2)\A N ) h , 

fc + V -(2)|B 4 ) fc 

■■ + V -(JV/2)|B JV ) fc . 



(C-21) 



(C-22) 



(C-23) 



(C-24) 



(C-25) 



(C-26) 



(C-27) 



(C-28) 



(C-29) 



Repeating the treatment similar to that presented in the 
odd-A^ case, we can reduce Gff, m hito Gfq tm which again 
obeys eq. (56). 
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